Modelling the impact of forest management and CO2-fertilisation on growth and demography in a Sitka spruce plantation

Afforestation and reforestation to meet ‘Net Zero’ emissions targets are considered a necessary policy by many countries. Their potential benefits are usually assessed through forest carbon and growth models. The implementation of vegetation demography gives scope to represent forest management and other size-dependent processes within land surface models (LSMs). In this paper, we evaluate the impact of including management within an LSM that represents demography, using both in-situ and reanalysis climate drivers at a mature, upland Sitka spruce plantation in Northumberland, UK. We compare historical simulations with fixed and variable CO2 concentrations, and with and without tree thinning implemented. Simulations are evaluated against the observed vegetation structure and carbon fluxes. Including thinning and the impact of increasing CO2 concentration (‘CO2 fertilisation’) gave more realistic estimates of stand-structure and physical characteristics. Historical CO2 fertilisation had a noticeable effect on the Gross Primary Productivity seasonal–diurnal cycle and contributed to approximately 7% higher stand biomass by 2018. The net effect of both processes resulted in a decrease of tree density and biomass, but an increase in tree height and leaf area index.


Supplementary Information: JULES-RED Model Description
The Joint UK Land Environment Simulator (JULES) Land Surface Model (LSM) is used to simulate the bio-geophysics of the land surface [1,2] across a landscape represented by tiles. JULES calculates the photosynthetic pathways for C3 and C4 using the Collatz et al., 1991[3] and , 1992 [4] models, respectively. In the JULES model, the soil hydrology is simulated using a discretisation of the Richards equation and the van Genuchten [5] water retention and hydraulic conductivity curves. As outlined in Clark, D. et al., 2011[2], tile carbon pools and litter fluxes were estimated for wood, leaf, and root, using Plant Functional Type (PFT) tile allometric and litter production rate equations. These relationships have been improved with evaluations against remotely sense and site comparisons, with more specific PFTs, such as Needleleaf Evergreen Tree (NET), being introduced [6,7]. JULES calculates the carbon and water fluxes at a half hourly timestep, with the vegetation dynamics (tree carbon, height and LAI), being updated daily.
The Robust Ecosystem Demography (RED) DGVM has been coupled into the JULES LSM [8]. JULES-RED partitions the number density, (kgC -1 m -2 ), of each Plant Functional Type (PFTs) into mass, (kgC), size classes and updates the size-structure by using equation (S1), a Fokker-Planck continuity equation of plant growth, (kgC yr -1 ), and mortality, (yr -1 ): Grid-box vegetation coverage, or tile fraction, is normally estimated by taking the integral of the product of number density and crown area, (m 2 ). However, this integral can potentially exceed the grid-box area when the number density is sufficiently high, such as in a Sitka spruce plantation [9]. Therefore, we have implemented a non-restricted "crown-area" fraction, CA , (equation (S2)) which is not truncated to 1: The "top-down" or grid-box vegetation fraction for PFT is given in equation (S3): where >1 is the canopy reduction from CA to be seen from above the canopy. The necessary reduction of >1, for each PFT, is estimated from the cumulative PFT sum in order of ascending height in the mass classes and PFT dimensions, until the total difference (∑ CA, − >1, ) across PFTs drops below 1. RED takes JULES inputs of carbon assimilate density ( ), which in JULES-RED is the difference of estimated NPP (Π NPP ) and the local litterfall (Λ LLF ) multiplied by the grid-box coverage (for PFT grid-box density), as demonstrated by equation (S4):

= (Π NPP − Λ LLF ). (S4)
The addition of canopy closure results in a simplistic method of limiting the overall growth rate as coverage is truncated to 1. Equation (S5) shows the recruitment dynamics. Recruitment assumes a fraction ( ) of the PFT assimilate is devoted to seedling reproduction: To simulate JULES-RED, we require radiative and meteorological forcings: air pressure, air temperature, specific humidity, precipitation, wind speed, downward short and long wave radiation.
To run with daily data JULES-RED also requires the daily temperature range. We also require soil property variables to run the van Genuchten model.
For simulating small initial trees an additional lower mass class of 0 = 0.1 kgC (normally 1.0 kgC) was added onto the lowest NET mass class, and we simulate 30 geometrically spaced mass classes up to 50,000 kgC (standard configuration maximum tree mass class in RED). Additionally, the allometric parameters for 0 = 0.23m 2 , ℎ 0 = 2.89m, bal,0 = 0.8 m 2 m −2 . Allometry follows a power-law relationship with mass (equation (S6)), with the power for height and balanced LAI ( bal ) being 0.25. Generally, these power-law allometries, especially for height, are less accurate for smaller than for larger trees [11]. For NET, we assume a small, reseed fraction of = 0.005, this is to mimic the low recruitment rate of trees in these dense planted stands of young Sitka spruce trees [14]. We assume a baseline mortality rate of = 0.01 yr −1 , this value is lower than the mean mortality rate for Sitka spruce presented in the Forestry Commission yield management booklet for an un-thinned stand of ages between 35-49 [15].